Introducción

Es fascinante poder conocer ciertos aspectos del mundo del software libre y de código abierto orientado al análisis de datos espaciales, en está oportunidad quiero mostrarte un pequeño entrenamiento integrando dos softwares super potentes dentro del campo de la geocomputación, el primero es SAGA, el cual es un software de Sistema de Información Geográfica de código abierto y el otro es R, el cual es un lenguaje de programación libre e interpretado y orientado a objetos, este mismo presenta un gran ecosistema de paquetes para una diversidad de ejes temáticos, entre ellas podemos mencionar a la estadística, la estadística espacial, la teledetección, web mapping, machine learning, entre otros.

Para este entrenamiento tenemos que tener ciertos requerimientos como tener instalado SAGA y R, de preferencia se recomienda usar la IDE de Rstudio para tener una forma más amigable al momento de programar con R.

Observación:

Instalación de los paquetes necesarios para nuestro entrenamiento

libs <- c('RSAGA', 'tidyverse','raster',
          'sp','sf', 'tmap', 'cptcity')
install.packages(libs, dependencies = TRUE)

Activación de los paquetes instalados

library(tidyverse) # Pkgs para ciencia de datos 
library(RSAGA)     # Pkgs para integrar SAGA con R
library(raster)    # Pkgs para manejar datos raster
library(raster)    # Pkgs para manejar datos vectoriales
library(tmap)      # Pkgs para elaborar mapas temáticos
library(cptcity)   # Pkgs para obtener paleta de colores

Indentificando el directorio de trabajo donde están los módulos de SAGA

Observación:

env <- rsaga.env(path = 'C:/Program Files/QGIS 3.10/apps/saga-ltr')
env <- rsaga.env()
## Search for SAGA command line program and modules... 
## Done
env
## $workspace
## [1] "."
## 
## $cmd
## [1] "saga_cmd"
## 
## $path
## [1] "/usr/bin"
## 
## $modules
## [1] "/usr/lib/x86_64-linux-gnu/saga"
## 
## $version
## [1] "2.3.1"
## 
## $cores
## [1] NA
## 
## $parallel
## [1] FALSE
## 
## $lib.prefix
## [1] "lib"

El proceso anterior nos ayuda a integrar SAGA con R, adicionalmente nos brinda una información acerca de la versión de SAGA con la cual estamos trabajando, la ruta de los módulos y alguna información adicional de interés.

Ahora vamos a proceder a visualizar las librerias disponibles de SAGA para trabajar dentro de R

libs <- rsaga.get.libraries(path = env$modules)
libs
##  [1] "climate_tools"           "contrib_perego"         
##  [3] "db_odbc"                 "db_pgsql"               
##  [5] "docs_html"               "docs_pdf"               
##  [7] "garden_3d_viewer"        "garden_fractals"        
##  [9] "garden_games"            "garden_learn_to_program"
## [11] "garden_webservices"      "grid_analysis"          
## [13] "grid_calculus_bsl"       "grid_calculus"          
## [15] "grid_filter"             "grid_gridding"          
## [17] "grid_spline"             "grid_tools"             
## [19] "grid_visualisation"      "imagery_classification" 
## [21] "imagery_maxent"          "imagery_opencv"         
## [23] "imagery_photogrammetry"  "imagery_segmentation"   
## [25] "imagery_svm"             "imagery_tools"          
## [27] "imagery_vigra"           "io_esri_e00"            
## [29] "io_gdal"                 "io_gps"                 
## [31] "io_grid_image"           "io_grid"                
## [33] "io_shapes_dxf"           "io_shapes"              
## [35] "io_table"                "io_virtual"             
## [37] "pj_georeference"         "pj_proj4"               
## [39] "pointcloud_tools"        "pointcloud_viewer"      
## [41] "shapes_grid"             "shapes_lines"           
## [43] "shapes_points"           "shapes_polygons"        
## [45] "shapes_tools"            "shapes_transect"        
## [47] "sim_cellular_automata"   "sim_ecosystems_hugget"  
## [49] "sim_erosion"             "sim_hydrology"          
## [51] "sim_ihacres"             "sim_qm_of_esp"          
## [53] "sim_rivflow"             "statistics_grid"        
## [55] "statistics_kriging"      "statistics_points"      
## [57] "statistics_regression"   "ta_channels"            
## [59] "ta_compound"             "ta_hydrology"           
## [61] "ta_lighting"             "ta_morphometry"         
## [63] "ta_preprocessor"         "ta_profiles"            
## [65] "ta_slope_stability"      "table_calculus"         
## [67] "table_tools"             "tin_tools"              
## [69] "tin_viewer"
# Usamos la función lenght para indentificar el total de liberías: 
length(libs)
## [1] 69

Ahora vamos extraer todos los módulos de las liberías usando un pequeño código gracias al paquete tidyverse

# Estos son los módulos para una libería específica ('ta_channels')
mod <- rsaga.get.modules(libs = 'ta_channels', env = env)

# Ahora vamos a crear una tabla con toda las liberías 
# y sus respectivos módulos:
rsaga_mod <- rsaga.get.modules(libs = libs,env = env) %>% 
  bind_rows(.id = 'librerias') %>% as_tibble()

# Solo visualizamos las 10 primeras filas de toda la tabla
head(rsaga_mod,10)
## # A tibble: 10 x 4
##    librerias      code name                                 interactive
##    <chr>         <dbl> <chr>                                <lgl>      
##  1 climate_tools     0 Multi Level to Surface Interpolation FALSE      
##  2 climate_tools     1 Multi Level to Points Interpolation  FALSE      
##  3 climate_tools     2 Earth's Orbital Parameters           FALSE      
##  4 climate_tools     3 Annual Course of Daily Insolation    FALSE      
##  5 climate_tools     4 Daily Insolation over Latitude       FALSE      
##  6 climate_tools     5 Monthly Global by Latitude           FALSE      
##  7 climate_tools     6 PET (after Hargreaves, Table)        FALSE      
##  8 climate_tools     7 Daily to Hourly PET                  FALSE      
##  9 climate_tools     8 PET (after Hargreaves, Grid)         FALSE      
## 10 climate_tools     9 Sunrise and Sunset                   FALSE

Ahora vamos a trabajar con el paquete raster para poder leer el Modelo de Elevación Digital (DEM) que está esta disponible en el siguiente enlace click aquí

dem <- raster('Insumos/DEM.tif')
plot(dem)

Es importante corregir algunos valores vacíos que presenta nuestro modelo de elevación digital, para esto vamos a usar el algoritmo de Fill sinks de Wang & Liu.

rsaga.fill.sinks(in.dem = 'Insumos/DEM.tif',
                 out.dem = 'Insumos/DEM_fill')
## Search for SAGA command line program and modules... 
## Done
## 
## 
## SAGA Version: 2.3.1
## 
## library path: /usr/lib/x86_64-linux-gnu/saga/
## library name: libta_preprocessor
## library     : Preprocessing
## tool        : Fill Sinks (Planchon/Darboux, 2001)
## author      : Copyrights (c) 2003 by Volker Wichmann
## processors  : 4 [4]
## 
## loading: DEM
## 
## 
## Parameters
## 
## 
## Grid system: 30; 558x 661y; 963390.8124x 8058312.3304y
## DEM: DEM
## Filled DEM: Filled DEM
## Minimum Slope [Degree]: 0.010000
## 
## Save grid: Insumos/DEM_fill.sgrd...
# Visualización de nuestro DEM corregido
tmap_mode('view')
raster('Insumos/DEM_fill.sdat') %>% 
  tm_shape() + 
  tm_basemap('Esri.WorldTopoMap') +
  tm_raster(style = 'cont',alpha = 0.75,
            palette = cpt(pal = 'ncl_topo_15lev'))

Ahora procedemos a calcular el mapa de pendientes, aspecto y de curvatura

rsaga.slope.asp.curv(in.dem = 'Insumos/DEM_fill.sdat',
                     out.slope = 'Salida/pendiente',
                     out.aspect = 'Salida/aspecto',
                     out.cgene = 'Salida/curvatura',
                     unit.slope = 1,
                     unit.aspect = 1,
                     method = 'maxslope')
## Search for SAGA command line program and modules... 
## Done
## 
## 
## SAGA Version: 2.3.1
## 
## library path: /usr/lib/x86_64-linux-gnu/saga/
## library name: libta_morphometry
## library     : Morphometry
## tool        : Slope, Aspect, Curvature
## author      : O.Conrad (c) 2001
## processors  : 4 [4]
## 
## loading: DEM_fill
## 
## 
## Parameters
## 
## 
## Grid system: 30; 558x 661y; 963390.8124x 8058312.3304y
## Elevation: DEM_fill
## Slope: Slope
## Aspect: Aspect
## General Curvature: General Curvature
## Profile Curvature: <not set>
## Plan Curvature: <not set>
## Flow Line Curvature: <not set>
## Method: maximum slope (Travis et al. 1975)
## Slope Units: degree
## Aspect Units: degree
## 
## Save grid: Salida/pendiente.sgrd...
## Save grid: Salida/aspecto.sgrd...
## Save grid: Salida/curvatura.sgrd...
# Ahora procedemos a visualizar el mapa de pendiente
raster('Salida/pendiente.sdat') %>% 
  tm_shape() + 
  tm_basemap('Esri.WorldTopoMap') +
  tm_raster(style = 'quantile',n = 7,palette = 'Spectral')
# Visualizando el mapa de aspectos
raster('Salida/aspecto.sdat') %>% 
  tm_shape() + 
  tm_basemap('Esri.WorldTopoMap') +
  tm_raster(style = 'cont',palette = cpt(pal = 'grass_gdd',rev = TRUE))
# Visualizando el mapa de curvatura general
raster('Salida/curvatura.sdat') %>% 
  tm_shape() + 
  tm_basemap('Esri.WorldTopoMap') +
  tm_raster(style = 'cont',palette = cpt(pal = 'grass_ramp'))

Ahora vamos a calcular el índice de humedad topográfica (TWI) y el índice de posición topográfica (TPI)

rsaga.get.usage(lib = 'ta_hydrology',module = 15)
## Search for SAGA command line program and modules... 
## Done
## library path: /usr/lib/x86_64-linux-gnu/saga/
## library name: libta_hydrology
## library     : Hydrology
## Usage: saga_cmd ta_hydrology 15 [-DEM <str>] [-WEIGHT <str>] [-AREA <str>] [-SLOPE <str>] [-AREA_MOD <str>] [-TWI <str>] [-SUCTION <double>] [-AREA_TYPE <str>] [-SLOPE_TYPE <str>] [-SLOPE_MIN <double>] [-SLOPE_OFF <double>] [-SLOPE_WEIGHT <double>]
##   -DEM:<str>             Elevation
##  Grid (input)
##   -WEIGHT:<str>          Weights
##  Grid (optional input)
##   -AREA:<str>            Catchment area
##  Grid (output)
##   -SLOPE:<str>           Catchment slope
##  Grid (output)
##   -AREA_MOD:<str>        Modified Catchment Area
##  Grid (output)
##   -TWI:<str>             Topographic Wetness Index
##  Grid (output)
##   -SUCTION:<double>      Suction
##  Floating point
##  Minimum: 0.000000
##  Default: 10.000000
##   -AREA_TYPE:<str>       Type of Area
##  Choice
##  Available Choices:
##  [0] absolute catchment area
##  [1] square root of catchment area
##  [2] specific catchment area
##  Default: 1
##   -SLOPE_TYPE:<str>      Type of Slope
##  Choice
##  Available Choices:
##  [0] local slope
##  [1] catchment slope
##  Default: 1
##   -SLOPE_MIN:<double>    Minimum Slope
##  Floating point
##  Minimum: 0.000000
##  Default: 0.000000
##   -SLOPE_OFF:<double>    Offset Slope
##  Floating point
##  Minimum: 0.000000
##  Default: 0.100000
##   -SLOPE_WEIGHT:<double> Slope Weighting
##  Floating point
##  Minimum: 0.000000
rsaga.geoprocessor(lib = 'ta_hydrology',module = 15,
                   param = list(DEM = 'Insumos/DEM_fill.sdat',
                                TWI = 'Salida/Twi' ),
                   env = env)
## 
## 
## SAGA Version: 2.3.1
## 
## library path: /usr/lib/x86_64-linux-gnu/saga/
## library name: libta_hydrology
## library     : Hydrology
## tool        : SAGA Wetness Index
## author      : (c) 2001 by J.Boehner, O.Conrad
## processors  : 4 [4]
## 
## loading: DEM_fill
## 
## 
## Parameters
## 
## 
## Grid system: 30; 558x 661y; 963390.8124x 8058312.3304y
## Elevation: DEM_fill
## Weights: <not set>
## Catchment area: Catchment area
## Catchment slope: Catchment slope
## Modified Catchment Area: Modified Catchment Area
## Topographic Wetness Index: Topographic Wetness Index
## Suction: 10.000000
## Type of Area: square root of catchment area
## Type of Slope: catchment slope
## Minimum Slope: 0.000000
## Offset Slope: 0.100000
## Slope Weighting: 1.000000
## 
## Create index: DEM_fill
## catchment area and slope...
## pass 1 (39594 > 0)
## pass 2 (9653 > 0)
## pass 3 (4059 > 0)
## pass 4 (1920 > 0)
## pass 5 (950 > 0)
## pass 6 (484 > 0)
## pass 7 (246 > 0)
## pass 8 (130 > 0)
## pass 9 (64 > 0)
## pass 10 (32 > 0)
## pass 11 (18 > 0)
## pass 12 (5 > 0)
## pass 13 (2 > 0)
## pass 14 (0 > 0)
## post-processing...
## topographic wetness index...
## Save grid: Salida/Twi...
rsaga.get.usage(lib = 'ta_morphometry',module = 18)
## Search for SAGA command line program and modules... 
## Done
## library path: /usr/lib/x86_64-linux-gnu/saga/
## library name: libta_morphometry
## library     : Morphometry
## Usage: saga_cmd ta_morphometry 18 [-DEM <str>] [-TPI <str>] [-STANDARD <str>] [-RADIUS_MIN <double>] [-RADIUS_MAX <double>] [-DW_WEIGHTING <str>] [-DW_IDW_POWER <double>] [-DW_IDW_OFFSET <str>] [-DW_BANDWIDTH <double>]
##   -DEM:<str>             Elevation
##  Grid (input)
##   -TPI:<str>             Topographic Position Index
##  Grid (output)
##   -STANDARD:<str>        Standardize
##  Boolean
##  Default: 0
##   -RADIUS_MIN:<double>   Radius
##  Value range
##   -RADIUS_MAX:<double>   Radius
##  Value range
##   -DW_WEIGHTING:<str>    Weighting Function
##  Choice
##  Available Choices:
##  [0] no distance weighting
##  [1] inverse distance to a power
##  [2] exponential
##  [3] gaussian weighting
##  Default: 0
##   -DW_IDW_POWER:<double> Inverse Distance Weighting Power
##  Floating point
##  Minimum: 0.000000
##  Default: 1.000000
##   -DW_IDW_OFFSET:<str>   Inverse Distance Offset
##  Boolean
##  Default: 1
##   -DW_BANDWIDTH:<double> Gaussian and Exponential Weighting Bandwidth
##  Floating point
##  Minimum: 0.000000
rsaga.geoprocessor(lib = 'ta_morphometry',module = 18,
                   param = list(DEM = 'Insumos/DEM_fill.sdat',
                                TPI = 'Salida/Tpi' ),
                   env = env)
## 
## 
## SAGA Version: 2.3.1
## 
## library path: /usr/lib/x86_64-linux-gnu/saga/
## library name: libta_morphometry
## library     : Morphometry
## tool        : Topographic Position Index (TPI)
## author      : O.Conrad (c) 2011
## processors  : 4 [4]
## 
## loading: DEM_fill
## 
## 
## Parameters
## 
## 
## Grid system: 30; 558x 661y; 963390.8124x 8058312.3304y
## Elevation: DEM_fill
## Topographic Position Index: Topographic Position Index
## Standardize: no
## Radius: 0.000000; 100.000000
## Weighting Function: no distance weighting
## 
## Save grid: Salida/Tpi...
# Visualizando el mapa de índice de humedad topográfica
raster('Salida/Twi.sdat') %>% 
  tm_shape() + 
  tm_basemap('Esri.WorldTopoMap') +
  tm_raster(style = 'cont',palette = cpt(pal = 'idv_blues'))
# Visualizando el mapa de índice de posición topográfica
raster('Salida/Tpi.sdat') %>% 
  tm_shape() + 
  tm_basemap('Esri.WorldTopoMap') +
  tm_raster(style = 'cont',palette = cpt(pal = 'neota_base_yellow'))

** Con RSAGA tambien podemos el orden de Strahler y la red de drenaje**

rsaga.get.usage(lib = 'ta_channels',module = 5)
## Search for SAGA command line program and modules... 
## Done
## library path: /usr/lib/x86_64-linux-gnu/saga/
## library name: libta_channels
## library     : Channels
## Usage: saga_cmd ta_channels 5 [-DEM <str>] [-DIRECTION <str>] [-CONNECTION <str>] [-ORDER <str>] [-BASIN <str>] [-SEGMENTS <str>] [-BASINS <str>] [-NODES <str>] [-THRESHOLD <num>]
##   -DEM:<str>         Elevation
##  Grid (input)
##   -DIRECTION:<str>   Flow Direction
##  Grid (optional output)
##   -CONNECTION:<str>  Flow Connectivity
##  Grid (optional output)
##   -ORDER:<str>       Strahler Order
##  Grid (optional output)
##   -BASIN:<str>       Drainage Basins
##  Grid (optional output)
##   -SEGMENTS:<str>    Channels
##  Shapes (output)
##   -BASINS:<str>      Drainage Basins
##  Shapes (output)
##   -NODES:<str>       Junctions
##  Shapes (optional output)
##   -THRESHOLD:<num>   Threshold
##  Integer
##  Minimum: 1
rsaga.geoprocessor(lib = 'ta_channels',module = 5,
                   param = list(DEM = 'Insumos/DEM_fill.sdat',
                                ORDER = 'Salida/Strahler_order',
                                SEGMENTS = 'Salida/red_drenaje'),
                   env = env)
## 
## 
## SAGA Version: 2.3.1
## 
## library path: /usr/lib/x86_64-linux-gnu/saga/
## library name: libta_channels
## library     : Channels
## tool        : Channel Network and Drainage Basins
## author      : O.Conrad (c) 2003
## processors  : 4 [4]
## 
## loading: DEM_fill
## 
## 
## Parameters
## 
## 
## Grid system: 30; 558x 661y; 963390.8124x 8058312.3304y
## Elevation: DEM_fill
## Flow Direction: <not set>
## Flow Connectivity: <not set>
## Strahler Order: Strahler Order
## Drainage Basins: <not set>
## Channels: Channels
## Drainage Basins: Drainage Basins
## Junctions: <not set>
## Threshold: 5
## 
## Flow Direction
## Stream Order
## Junctions
## Drainage Basins
## Vectorising Grid Classes
## 
## 
## Parameters
## 
## 
## Grid system: 30; 558x 661y; 963390.8124x 8058312.3304y
## Grid: Drainage Basins
## Polygons: Polygons
## Class Selection: all classes
## Vectorised class as...: one single (multi-)polygon object
## Keep Vertices on Straight Lines: no
## 
## class identification
## Create index: Drainage Basins
## edge detection
## edge collection
## Channels
## Save grid: Salida/Strahler_order...
## Save shapes: Salida/red_drenaje...

Finalmente visualizamos nuestros resultados

** Strahler order

raster('Salida/Strahler_order.sdat') %>% 
  tm_shape() + 
  tm_basemap('Esri.WorldTopoMap') +
  tm_raster(style = 'fixed',breaks = seq(0,5,1),palette = '-viridis')

Red de drenaje

read_sf('Salida/red_drenaje.shp') %>% 
  tm_shape() + 
  tm_basemap('Esri.WorldTopoMap') +
  tm_lines()